Johnson Hsieh
April, 2015 @ DSP A1 coures
Contact me| Data Analyst Training Program @ DSP
ubike <- fread("data/ubikeweatherutf8.csv",
showProgress = interactive(), data.table = FALSE)
ubike1 <- filter(ubike, sno==1) %>%
mutate(sbi.range=max.sbi-min.sbi) %>%
mutate(is.rushhours=cut(hour, breaks=c(0, 8, 10, 17, 20, 24),
labels = c(0,1,0,1,0), right=FALSE)) %>%
mutate(is.weekday=ifelse(strftime(date, "%u") < 6, 1, 0))
tab1 <- filter(ubike1, is.rushhours==1, is.weekday==1) %>%
group_by(tot) %>%
summarise(min(sbi.range), mean(sbi.range), max(sbi.range))
df1 <- group_by(ubike1, date, hour) %>%
summarise(rate.sbi=mean(avg.sbi)/tot) %>%
group_by(hour) %>%
summarise(rate.sbi=mean(rate.sbi))
thm <- function() {
theme_gray(base_family = "STHeiti") +
theme(text=element_text(size=18))
}
ggplot(df1, aes(x=hour, y=rate.sbi)) +
geom_bar(stat="identity") +
ggtitle("捷運市政府站") +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank()) +
scale_y_continuous(labels = percent)
rate.sbi) 為 平均車輛數 / 總停車格數 (avg.sbi/tot)is.rain),當該時段累積雨量大於1mm訂為雨天,反之為晴天 依 日期 (date)、時間 (hour)、是否下雨 (is.rain) 做分組 (group_by) 計算 平均有車率 (rate.sbi = mean(avg.sbi/tot)),得到下表:df2 <- filter(ubike, sno==1) %>%
mutate(is.rain=rainfall>1) %>%
mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE),
labels = c("晴天","雨天"))) %>%
select(date, hour, tot, avg.sbi, avg.bemp, temp, is.rain) %>%
group_by(date, hour, is.rain) %>%
summarise(rate.sbi=mean(avg.sbi)/tot) %>%
group_by(hour, is.rain) %>%
summarise(rate.sbi=mean(rate.sbi))
kable(df2[17:22,], format="html", digits=3, align="c")
| hour | is.rain | rate.sbi |
|---|---|---|
| 8 | 晴天 | 0.431 |
| 8 | 雨天 | 0.658 |
| 9 | 晴天 | 0.288 |
| 9 | 雨天 | 0.576 |
| 10 | 晴天 | 0.254 |
| 10 | 雨天 | 0.547 |
首先用長條圖 (bar chart) 來探索這份報表,當欄位大於二時,將依賴顏色做區隔,一般而言長條圖有以下變型: - Stack plot (堆疊圖) - Dodge plot () - Fill plot (相對堆疊) - Facet panels - Pyramid (金字塔圖)
geom_bar(stat="identity", position="stack")ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
geom_bar(stat="identity") +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank())
geom_bar(stat="identity", position="dodge")ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
geom_bar(stat="identity", position="dodge") +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank())
geom_bar(stat="identity", position="fill")ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
geom_bar(stat="identity", position="fill") +
labs(x="時間", y="相對有車率") +
thm() +
theme(legend.title=element_blank())
facet_grid(y~.) or facet_grid(.~x)ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
geom_bar(stat="identity", position="dodge") +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank()) +
facet_grid(is.rain~.)
ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
geom_bar(stat="identity", position="dodge") +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank()) +
facet_grid(.~is.rain)
filter(df2, is.rain=="晴天"), and coord_flip()ggplot(df2, aes(x=hour,y=rate.sbi, fill=is.rain)) +
geom_bar(data=filter(df2, is.rain=="晴天"), stat="identity") +
geom_bar(aes(y=rate.sbi*(-1)), data=filter(df2, is.rain=="雨天"),
stat="identity") +
scale_y_continuous(breaks=seq(from=-1, to=1, by=0.1),
labels=abs(seq(-1, 1, 0.1))) +
labs(x="時間", y="有車率") +
theme(legend.title=element_blank()) +
coord_flip() + thm()
由於x軸的單位為時間,所以可以改用折線圖 (line chart) 來做探索,本例就比較晴、雨天有車率的差異而言,line chart 比 bar chart 更具視覺上優勢。
ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
geom_line(size=1) +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank())
stat_smooth()ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
#geom_line() +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank()) +
stat_smooth(size=1)
coord_polar()ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
#geom_line() +
labs(x="時間", y="有車率") +
thm() +
theme(legend.title=element_blank()) +
stat_smooth(size=1) + coord_polar()
熱點圖 (heatmap) 是用顏色深淺呈現數值大小的視覺化。
Hint: geom_tile()
ggplot(df2, aes(x=hour, y=is.rain, fill=rate.sbi)) +
geom_tile() +
scale_fill_gradient(name="有車率", low="white", high="midnightblue") +
labs(x="時間", y="天氣") +
thm()
平行座標圖 (Parallel coordinate plot) 多用於呈現多欄位的資料視覺化,強調欄位的順序性,特別適合用在因果關係的陳述。譬如:行業別 -> 是否上DSP課程 -> 職場表現。
Hint: library(GGally) and ggparcoord()
df2 <- mutate(df2, rain=as.numeric(is.rain)-1)
ggparcoord(as.data.frame(df2), columns = c(1,4,3), groupColumn = 2) +
thm() + theme(legend.title=element_blank())
ggparcoord(data = iris, columns = 1:4, groupColumn = 5,
title = "Parallel Coordinate Plot for the Iris Data") + thm()
rate.sbi) 為 平均車輛數 / 總停車格數 (avg.sbi)/totrate.used) 為 (最大車輛數 - 最小車輛數) / 總停車格數 (max.sbi-min.sbi)/tot
is.weekday)、是否下雨 (is.rain)、時間 (hour)、有車率 (rate.sbi)、使用率 (rate.used) 的關係dplyr中filter, mutate, group_by, summarise等函數將資料整理成下表:df3 <- filter(ubike1, sno==1) %>%
mutate(is.rain=rainfall>1) %>%
mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE),
labels = c("晴天","雨天"))) %>%
mutate(is.weekday=strftime(date, "%u")<6) %>%
mutate(is.weekday=factor(is.weekday, levels=c(FALSE, TRUE),
labels=c("平日","假日"))) %>%
mutate(is.rushhours=cut(hour, breaks=c(0, 4, 7, 24), right=FALSE)) %>%
group_by(date, is.weekday, is.rushhours, hour, is.rain) %>%
summarise(rate.sbi=mean(avg.sbi)/tot, rate.used=mean(max.sbi-min.sbi)/tot)
kable(head(df3), format="html", digits=3, align="c")
| date | is.weekday | is.rushhours | hour | is.rain | rate.sbi | rate.used |
|---|---|---|---|---|---|---|
| 2014-12-08 | 假日 | [7,24) | 15 | 雨天 | 0.536 | 0.033 |
| 2014-12-08 | 假日 | [7,24) | 16 | 雨天 | 0.545 | 0.044 |
| 2014-12-08 | 假日 | [7,24) | 17 | 雨天 | 0.547 | 0.039 |
| 2014-12-08 | 假日 | [7,24) | 18 | 雨天 | 0.573 | 0.044 |
| 2014-12-08 | 假日 | [7,24) | 19 | 雨天 | 0.556 | 0.100 |
| 2014-12-08 | 假日 | [7,24) | 20 | 雨天 | 0.737 | 0.100 |
散佈圖 (scatterplot) 是比較兩數值變數最直覺的視覺化,主要是觀察兩變數間是否存在特殊的趨勢、群聚現象、離群值。以下方式常用來輔助散佈圖的探索:
ggplot(df3, aes(x=rate.used, y=rate.sbi)) + geom_point() +
labs(x="使用率", y="有車率") + thm()
ggplot(df3, aes(x=rate.used, y=rate.sbi)) + geom_point() +
labs(x="使用率", y="有車率") + thm()+
geom_vline(xintercept=0.4, lty=2) +
geom_hline(yintercept=0.4, lty=2)
aes(colour=is.rushhours)ggplot(df3, aes(x=rate.used, y=rate.sbi)) +
labs(x="使用率", y="有車率") + thm() +
geom_point(aes(colour=is.rushhours), position="jitter")
aes(colour=is.rushhours, shape=is.weekday)ggplot(df3, aes(x=rate.used, y=rate.sbi)) +
geom_point(aes(colour=is.rushhours, shape=is.weekday), position="jitter") +
facet_grid(is.weekday~is.rushhours) +
labs(x="使用率", y="有車率") +
thm() +
theme(legend.title=element_blank())
df4 <- filter(df3, hour>6 & hour < 10)
ggplot(df4, aes(x=rate.used, y=rate.sbi)) +
geom_point(aes(colour=paste(as.character(is.weekday),
as.character(is.rain), sep="-"))) +
ggtitle("每日7-9點 YouBike使用狀況") +
labs(x="使用率", y="有車率") +
facet_grid(is.rain~is.weekday) +
thm() +
theme(legend.title=element_blank())
tmp <- group_by(ubike, sno, sna, sarea, lat, lng) %>% distinct
dist <- round(distm(x=tmp[, c("lng","lat")])[,1])
df5 <- tmp %>% select(sno, sna, sarea, lat, lng) %>%
cbind(dist) %>% arrange(dist) %>% top_n(10, wt = -dist)
研究完單一場站之後,試著探索市府站與鄰近場站的關係,此時需要透過經緯度計算場站與場站之間的距離。透過geosphere套件中的distm函數可以批次計算所有場站之間的兩兩距離,整理得到下表,離捷運市府站最近的場站依序是台北市政府 (438m), 興雅國中 (484m)…。
Hint: library(geosphere), distm, group_by, distinct
kable(df5, format="html", digits=3, align="c")
| sno | sna | sarea | lat | lng | dist |
|---|---|---|---|---|---|
| 1 | 捷運市政府站(3號出口) | 信義區 | 25.041 | 121.568 | 0 |
| 3 | 台北市政府 | 信義區 | 25.038 | 121.565 | 438 |
| 5 | 興雅國中 | 信義區 | 25.037 | 121.569 | 484 |
| 25 | 永吉松信路口 | 信義區 | 25.045 | 121.572 | 659 |
| 6 | 世貿二館 | 信義區 | 25.035 | 121.566 | 718 |
| 150 | 松德公園 | 信義區 | 25.037 | 121.573 | 734 |
| 138 | 捷運永春站(2號出口) | 信義區 | 25.041 | 121.575 | 754 |
| 8 | 世貿三館 | 信義區 | 25.035 | 121.564 | 759 |
| 113 | 仁愛逸仙路口 | 信義區 | 25.038 | 121.561 | 763 |
| 4 | 市民廣場 | 信義區 | 25.036 | 121.562 | 778 |
利用ggmap套件導入google map作為底圖將場站位置標示出來。
Hint: library(ggmap), map <- get_map("Taipei"); ggmap(map), geom_point
df5$is.cityhall <- factor(c(1, rep(0, 9)), levels=1:0)
map <- get_map(location=c(lon=df5$lng[1], lat=df5$lat[1]) , zoom = 15)
ggmap(map) + thm() +
geom_point(data=df5, aes(x=lng, y=lat, colour=is.cityhall), size=5) +
geom_text(data=df5, aes(x=lng, y=lat, label=sna, colour=is.cityhall),
position="jitter", vjust=-1, hjust=0.5, size=4, family="STHeiti") +
theme(legend.position="none") + scale_color_brewer(palette="Set1")
練習用geom_point(size=tot)來改變場站標示的大小。
df5 <- group_by(tmp) %>% select(sno, tot) %>%
right_join(df5, by="sno") %>%
`[`(c(1, 3, 4, 2, 5, 6, 7, 8))
map <- get_map(location=c(lon=df5$lng[1], lat=df5$lat[1]) ,
maptype = "roadmap", zoom = 15)
ggmap(map) + thm() +
geom_point(data=df5, aes(x=lng, y=lat, colour=is.cityhall, size=tot)) +
theme(legend.position="none") + scale_color_brewer(palette="Set1") +
scale_size(range = c(3,12))
觀察鄰近捷運市府站的10個YouBike場站,每一天 有車率 與 使用率的狀況。以有車率為例,透過觀察可以發現{興雅國中, 台北市政府, 市民廣場, 世貿三館, 世貿二館} 時間分佈有相似的狀況,{永吉松信路口, 松德公園, 捷運永春站} 也有相似的情況,而捷運市府站介於兩群之間,仁愛逸仙路口則是一枝獨秀。
tmp1 <- filter(ubike, sno%in%df5$sno) %>%
mutate(is.rain=rainfall>1) %>%
mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE),
labels = c("晴天","雨天"))) %>%
mutate(is.weekday=strftime(date, "%u")<6) %>%
mutate(is.weekday=factor(is.weekday, levels=c(FALSE, TRUE),
labels=c("平日","假日"))) %>%
mutate(is.rushhours=cut(hour, breaks=c(0, 4, 7, 24), right=FALSE)) %>%
group_by(date, sno, sna, is.weekday, is.rushhours, is.rain, hour, tot) %>%
summarise(rate.sbi=mean(avg.sbi)/tot, rate.used=mean(max.sbi-min.sbi)/tot)
df6 <- tmp1 %>%
filter(is.weekday=="平日", is.rain=="晴天") %>%
group_by(sno, sna, sna, hour) %>%
summarise(rate.sbi=mean(rate.sbi), rate.used=mean(rate.used))
ggplot(df6, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() +
theme(legend.position="bottom") +
scale_fill_gradient(name="有車率", low="white", high="lawngreen") +
labs(x="時間", y="") +
theme(axis.text = element_text(size = 13, color="darkgreen"))
ggplot(df6, aes(x=hour, y=sna, fill=rate.used)) + geom_tile() + thm() +
theme(legend.position="bottom") +
scale_fill_gradient(name="使用率", low="white", high="Navy") +
labs(x="時間", y="") +
theme(axis.text = element_text(size = 13, color="darkblue"))
當heatmap的x軸或y軸為類別變數時,可以經由階層分群法 (hierarchical clustering) 做行或列的排序。
dcast)hclust)ggdendrogram)取得排序 (order)
Hint: library(reshape2), library(ggdendro)
dat <- dcast(df6, sna~hour, value.var="rate.sbi")
rownames(dat) <- dat[,1]
dat <- dat[,-1]
kable(dat[,8:13], format = "html", digits = 3, row.names = TRUE)
| 7 | 8 | 9 | 10 | 11 | 12 | |
|---|---|---|---|---|---|---|
| 捷運市政府站(3號出口) | 0.589 | 0.553 | 0.504 | 0.463 | 0.479 | 0.484 |
| 捷運永春站(2號出口) | 0.276 | 0.229 | 0.188 | 0.192 | 0.160 | 0.156 |
| 仁愛逸仙路口 | 0.228 | 0.184 | 0.157 | 0.151 | 0.137 | 0.192 |
| 世貿二館 | 0.148 | 0.216 | 0.372 | 0.522 | 0.760 | 0.767 |
| 世貿三館 | 0.146 | 0.197 | 0.270 | 0.397 | 0.469 | 0.465 |
| 市民廣場 | 0.204 | 0.257 | 0.316 | 0.389 | 0.487 | 0.510 |
| 松德公園 | 0.267 | 0.218 | 0.156 | 0.095 | 0.059 | 0.039 |
| 台北市政府 | 0.251 | 0.242 | 0.297 | 0.386 | 0.497 | 0.625 |
| 興雅國中 | 0.192 | 0.287 | 0.418 | 0.543 | 0.627 | 0.686 |
| 永吉松信路口 | 0.280 | 0.245 | 0.166 | 0.092 | 0.058 | 0.049 |
hc.sna <- hclust(dist(dat))
ggdendrogram(hc.sna, rotate = TRUE) + thm() + labs(x="", y="")
# hc.sna$order
sna.order <- data.frame(order=1:10, sna=hc.sna$labels[hc.sna$order])
kable(sna.order, format = "html")
| order | sna |
|---|---|
| 1 | 仁愛逸仙路口 |
| 2 | 捷運永春站(2號出口) |
| 3 | 松德公園 |
| 4 | 永吉松信路口 |
| 5 | 捷運市政府站(3號出口) |
| 6 | 世貿三館 |
| 7 | 市民廣場 |
| 8 | 世貿二館 |
| 9 | 台北市政府 |
| 10 | 興雅國中 |
df7 <- df6
df7$sna <- factor(df7$sna, levels=(sna.order[,2]))
ggplot(df7, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() +
theme(legend.position="bottom") +
scale_fill_gradient(name="有車率", low="white", high="lawngreen") +
labs(x="時間", y="") +
theme(axis.text = element_text(size = 13, color="darkgreen"))
hc.hour <- hclust(dist(t(dat)))
ggdendrogram(hc.hour) + thm() + labs(x="", y="")
hour.order <- data.frame(order=1:24, sna=hc.hour$labels[hc.hour$order])
# kable(hour.order, format = "html")
df7$hour <- factor(df7$hour, levels=(hour.order[,2]))
ggplot(df7, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() +
theme(legend.position="bottom") +
scale_fill_gradient(name="有車率", low="white", high="lawngreen") +
labs(x="時間", y="") +
theme(axis.text = element_text(size = 13, color="darkgreen"))
dat <- dcast(df6, sna~hour, value.var="rate.used")
rownames(dat) <- dat[,1]
dat <- dat[,-1]
hc.sna <- hclust(dist(dat))
hc.hour <- hclust(dist(t(dat)))
df8 <- df6
df8$sna <- factor(df8$sna, levels = hc.sna$labels[hc.sna$order])
df8$hour <- factor(df8$hour, levels = hc.hour$labels[hc.hour$order])
ggplot(df8, aes(x=hour, y=sna, fill=rate.used)) + geom_tile() + thm() +
theme(legend.position="bottom") +
scale_fill_gradient(name="使用率", low="white", high="Navy") +
labs(x="時間", y="") +
theme(axis.text = element_text(size = 13, color="darkblue"))
平行座標圖常用來展示不同群組在諸多變數間的差異性,當群組分類方式未知時,可以利用機器學習 (machine learning) 中的非監督式學習 (unsupervised learning),幫資料做分群。分群之後再藉由平行座標圖來呈現資料的脈絡。
tot)、有車率 (rate.sbi)、使用率 (rate.used) 三個變數做分群
tmp2 <- filter(tmp1, is.weekday=="平日", is.rain=="晴天", hour>6 & hour<22) %>%
group_by(sno, sna, tot) %>%
summarise(rate.sbi=mean(rate.sbi), rate.used=mean(rate.used))
km <- kmeans(tmp2[,3:5], 3)
km
K-means clustering with 3 clusters of sizes 1, 4, 5
Cluster means:
tot rate.sbi rate.used
1 180.0 0.3260354 0.1156761
2 65.0 0.3684294 0.1783035
3 35.2 0.2290803 0.2471085
Clustering vector:
[1] 1 3 2 2 2 2 3 3 3 3
Within cluster sum of squares by cluster:
[1] 0.0000 300.0137 92.8441
(between_SS / total_SS = 97.8 %)
Available components:
[1] "cluster" "centers" "totss" "withinss"
[5] "tot.withinss" "betweenss" "size" "iter"
[9] "ifault"
df9 <- group_by(tmp2) %>%
transmute(sna, tot, rate.sbi, rate.used,
group=factor(km$cluster)) %>%
arrange(group)
ggparcoord(as.data.frame(df9), columns = c(1,2,3,4), groupColumn = 5,
scale="uniminmax") +
geom_line(size=1) + thm() + theme(legend.title=element_blank()) +
scale_x_discrete(labels=c("場站","總停車格","有車率","使用率")) +
labs(x="", y="")
kable(df9, format="html", digits=3, align="c")
| sna | tot | rate.sbi | rate.used | group |
|---|---|---|---|---|
| 捷運市政府站(3號出口) | 180 | 0.326 | 0.116 | 1 |
| 市民廣場 | 60 | 0.367 | 0.137 | 2 |
| 興雅國中 | 60 | 0.384 | 0.163 | 2 |
| 世貿二館 | 80 | 0.410 | 0.257 | 2 |
| 世貿三館 | 60 | 0.313 | 0.156 | 2 |
| 台北市政府 | 40 | 0.361 | 0.243 | 3 |
| 永吉松信路口 | 30 | 0.200 | 0.266 | 3 |
| 仁愛逸仙路口 | 38 | 0.180 | 0.192 | 3 |
| 捷運永春站(2號出口) | 30 | 0.234 | 0.356 | 3 |
| 松德公園 | 38 | 0.170 | 0.179 | 3 |